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in the wave functions of the colliding nuclei. We argue that the dominant production mechanism in proton-nucleus 

O^l and nucleus-nucleus collisions for heavy nuclei is different from the one in hadron-hadron interactions. We find 
that the rapidity distribution of primary J/ip production is more peaked around midrapidity than the analogous 
distribution in elementary pp collisions. We discuss the consequences of this fact on the experimentally observed J/ip 

0^ suppression in Au — Au collisions at RHIC energies. 
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1. Introduction 

Understanding the mechanism of J/ip production has been a challenge for over three decades. Despite a relatively large 
charm quark mass, the binding energy of J/ip is quite small; therefore non-pcrturbative corrections can be important. 
This is a likely cause of the difficulties encountered by perturbative QCD in describing the differential J/ip production 
cross section and its polarization. A significant effort has been invested into attempts to uncover the mysteries of 
J/ip production. Still, when confronted with the experimental data the existing approaches encounter problems that 
have to be cured by the introduction of additional adjustable parameters encoding the poorly understood dynamics 
(for a recent review see [1]). 

In this paper we develop a new approach to the J/ip production in nuclear reaction suggested by two of us in 
[2]. It was argued in Rcf. [2] that at high energies the dynamics of J/ip production is determined mostly by the 
strength of the coherent quasi-classical fields of the nucleus. This approach yielded a reasonable description of the 
experimental data, and in particular shone some light on the possible origin of xp scaling in p(d)A collisions observed 
in the data from CERN[3], FNAL[4] and RHIC[5]. 

In the context of high energy nuclear physics, it is important to understand well the mechanism of J/ip production 
since J/ip suppression in heavy ion collisions could serve as a signal of the Quark-Gluon Plasma [6]. Motivated by 
the urgent necessity to understand the cold nuclear effects on the J/ip production in heavy ion collisions, we shall 
calculate the inclusive J/ip production in heavy ion collisions. In this paper we take into account only the cold 
nuclear matter effects neglecting completely the dynamical effects leading to the possible formation of the Quark 
Gluon Plasma. 
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A systematic approach to the particle production in heavy ion collisions at high energies has been developed 
in Refs. [7, 8, 9, 10] and will be referred to as the KLN approach. KLN assumes that the wave functions of the 
colliding nuclei can be described as the Color Glass Condensate (CGC) [12, 13] for which the most characteristic 
property is the saturation of the parton density [14, 15, 12, 16]. These ideas have passed the first check against 
the RHIC experimental data on multiplicities and rapidity distributions (see Refs. [7, 8, 9, 10, 11, 17]) and, in this 
paper, we confront them with the rapidity distribution of J/ip mesons at RHIC. Our study in this paper is based on 
the following observations made in Rcf. [2]: (1) The mechanism of J/ip production in hadron- nucleus collisions is 
different from the one in the hadron-hadron interactions; (2) Inclusive J /ip production and inclusive cc production 
are dominated by different distance scales. 

To explain our main idea, consider the J ftp production in hadron-hadron collisions. The leading contribution is 
given by the two-gluon fusion 

G + G -> J/ip + gluon. (1.1) 

This process is of the order O(a^): the partonic sub-process is of the order a 3 ; two additional powers of a s arise 
from attaching the initial gluons to the colliding hadrons. The thrcc-gluon fusion 

G + G + G -> J/ip (1.2) 

is parametrically suppressed as it is proportional to 0(a 6 s ). However, in hadron-nucleus collisions two of the initial 
state gluons can be attached to the nucleus. This brings in an additional enhancement by A 1 / 3 . Since in the quasi- 
classical approximation a^A 1 ^ ~ 1 we find that the thrcc-gluon fusion of (1.2) is actually enhanced by l/a s as 
compared to (1.1). Similar conclusion holds for heavy ion collisions. 

A particularly helpful insight into the nature of the contribution (1.2) is obtained if we note that three-gluon 
contribution (1.2) is suppressed as compared to the two-gluon one (1.1) by an additional factor r 2 , where (2to c ) _1 < 
r < (2m c a s )~ 1 . This factor arises since we need to have three gluons in the area of the order of r 2 . In other words, it 
means that this reaction originates from the higher-twist contribution. However, in the hadron - nucleus interactions 
the higher-twist contribution appears always in the dimcnsionless combination r 2 Q 2 with the saturation scale Q s . 
The saturation scale is proportional to A 1 / 3 which compensates for the smallncss of r. The dominance of the higher 
twist process (1.2) is main idea of [2] and we are going to develop it in this paper in the case of heavy ion collisions. 
Various aspects of multi-parton interactions generating the higher twist effects in J/ip production were considered 
previously in [18, 19, 20, 22, 23, 24, 21, 26, 27, 28]. 

The paper is structured as follows: in Section 2 we describe the production of cc pairs in proton- nucleus and 
nucleus- nucleus collisions; in Sections 3 and 4 we calculate the inclusive J/ip production cross-section in these 
collisions. We present our numerical results in Section 5. 



2. Warm-up: inclusive production of cc pair with fixed relative momentum 
2.1 Hadron-hadron collisions 

The process of cc production in a hadron-hadron collision is shown in Fig. 1. The corresponding cross section 
integrated over the transverse momentum of the pair, denoted by p = p (p 2 = p 2 ) in Fig. 1, is equal to the square of 
the diagram in this figure and can be written as follows 

dcrjpp) 1 1 47T 2 q s f d 2 h 2 f d 2 l 2 , 2 . 

dYd^k = (2^2(^-1) N c M —M*i,h) J ^o(^k) x 

d 2 rdz^ G (h,r,z)(l-e u -^) e - **** f d 2 r' V G (h,r' , z) (l - e'^A e l ^ r - (2.1) 
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Figure 1: The process of inclusive cc production with fixed relative momentum (k) in hadron-hadron collision. 

where x\ — {m c ,t + ^c,t) eY /V s an d £2 = {jnc,t + m^) eT Y j \fs , with 1" being the rapidity of quark and antiquark 
pair in the center-of-mass frame, m 2 t — m 2 + k 2 , m| t = m 2 + fc?,, s, s' are the quark and anti-quark helicities and 
A is the gluon polarization. The function 0g(:e,Z 2 ) is the probability to find a gluon with given x and transverse 
momentum I. It is related to the gluon distribution function xG(x,Q 2 ) as 

Q 2 

xG(x,Q 2 ) = J dl 2 cj)(x,l 2 ) (2.2) 

The factor 2 in front of (2.1) is a consequence of the s-channcl unitarity by which the inelastic cross section equals 
twice the imaginary part of the elastic scattering amplitude. Eq. (1) is written in the /c^-factorization approach 
which is believed to be valid in hadron-hadron collisions at not too high energies [29, 30, 31]. 
It is convenient to introduce the cross section for dipole-hadron interaction in the form [32] 

= ^ / ^(l- e ^i)0(,,I 2 ). (2.3) 

In the DGLAP approximation the dominant contribution to the integral over I comes from the region Ir < 2 where it 
picks up the leading logarithmic contribution. Integrating first over all directions of the vector I and then expanding 
the resulting Bessel function yields: 

. 87r 2 cv s f 2 ^ r dl 1 9,9 ,9, a> s ir 2 9 „, . , o> ,„ , N 

r 2 ) « J j -r 2 l 2 cj>(x, I 2 ) = r 2 xG(x, 4/r 2 ) , (2.4) 

where we used (2.2). 

Using (2.2) and (2.3) it is easy to rewrite (2.1) in the following form 

da(pp) 1 1 ^ 2 

d 2 r * G (m c , r, z = 1/2) e~ l J** f d 2 r' ** G (m c , r 1 , z = 1/2) e l a in (x 2 , r, r') (2.5) 
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where 

Vin{x2,r,r') = a{x 2 ,r 2 ) + a(x 2 ,r' 2 ) - a(x 2 , (r - r') 2 ) . (2.6) 

(The a notation is used to distinguish the dipole cross section defined in (2.6) from the inclusive heavy quark- 
antiquark inelastic cross section we discuss later, see (2.12)). In derivation of (2.5) we took into account only the 
DGLAP contribution to x%G (xi, m 2 .) and we treated the c-quark as a non-relativistic particle with z — 1/2. All 
these simplifications are not important for our main results but allow for a more compact notations. 

The gluon light-cone wave function is well-known (Refs. [33, 34, 35]). It has the simplest form for z = 1/2 and 
ii/4 < ra 2 , namely (see (A. 26)) 



* G (m c ,r,z= 1/2) = iL 

Z7T 



i m c K\ (rm c ) A 5 S . S > + Ko(rm c ) s m c (l + s\)5 s .- s > 



r 



(2.7) 



$ G (m c ,r,r',z = 1/2) = — ^ 2(N l_ 1) Yl *g(™ c , r, z = l/2)tf£(m c , r', z = 1/2) 



a s m 2 c 



(2tt) J 



(2tt)3 2(N 2 - 1) a 

\ f ' T 

■x^f-Ki (rm c ) Ki (r m c ) + K (rm c ) K (r'm c ) 

2 TV 



(2. 



where t a is the Gell-Mann matrix and e is the polarization vector. With these definition we can write (2.5) as 



Xl G{ Xl , m 2 c ) / d 2 r / d 2 r' $ G (m c , r, /, z = 1/2) e l s(f -r)-* r? r >) ( 2 . 9 ) 



da{pp) 
dYd 2 k 

In Appendix we give a detailed derivation of these formulas. 
2.2 Hadron heavy nucleus collisions 

Production of quark-ant iquark pairs in high energy proton-nucleus collisions and in DIS both in the quasi-classical 
approximation of McLerran-Venugopalan model [12] (summing powers of a^.A 1 ^ 3 ) and including quantum small- a; 
evolution (summing powers of a s In i) has been calculated in Ref. [36, 37]. This process has been also considered 
by other authors [38, 39, 40] who obtained similar, though less general, results. Phcnomcnological applications have 
been addressed in details in [41, 42]. Using the results of [36, 37, 42] it is not difficult to generalize the formulae of 
the previous subsection for the case of pA collisions. The details arc given in Appendix. Here we present a derivation 
that emphasizes the key physical issues. 

As one can see in Fig. 2 the quark-antiquark pair production in hadron-nucleus interaction includes an additional 
clastic scattering of dipolcs with sizes r and r' as well as inelastic interaction at points Zi, which are the longitudinal 
coordinates of nuclcons in the nucleus 1 . To include both processes we need to modify (2.1) in the following way 



dainipA) 
dYd 2 kd 2 b 



Xl G{xi,m 2 c ) / d 2 r / d 2 r' $ G (m c ,r,r', z = 1/2) e^t'-rM 



2R A 





oo 



X 



pa m (x 2 ,r, r') dz e -[^2,r 2 )+a(x 2 ,r' 2 )] P 2 r a 

dz nP n al l {x^r,r') (2.10) 





p2R A p 


/ dzi .. 


. / rfz„_i / 


'ZQ 


J Zn—2 ** Z 



where p is the density of the nucleons in a nucleus and Ra is the nucleus radius. For brevity we wrote (2.10) for a 
cylindrical nucleus. In Sec. 5 we perform numerical analyses with realistic nuclear density distributions. 



1 Note that z (z 1 ) appearing in (2.7), (2.8) etc. denote the fraction of the gluon's light-cone momentum carried by the c-quark in the 
(complex conjugated) amplitude. Zj's with i = 0, 1, 2, ... in Fig. 2 etc. denote the longitudinal coordinates of nuclcons in the nucleus. 
These are two completely unrelated variables. 
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Figure 2: The process of inclusive cc production with fixed relative momentum (k) in hadron-nucleus collision. 



The factor exp { — [cr(x2, r 2 ) + a(x2 7 r' 2 )]p2R A } in (2.10) describes the fact that neither dipole with the size r 
(in the amplitude) nor dipole with the size r' (in the complex conjugate amplitude) interacts inclastically with the 
nucleons of the nucleus between the points and Z\ as well as between any other pair of points Zi and Zj_ i, where 
i = 0, 1, . . . ,n — 1. In deriving Eq. (2.10) we assumed that the collision energy is high enough so that the quark- 
antiquark pair is produced long before it starts to interact with the nucleus. This corresponds to coherent interaction 
of the pair with all nucleons of the target nucleus. It has been demonstrated in [41, 42] that this is indeed the case 
for charm quark production in central and forward 2 rapidities at RHIC and LHC. 

We assume in (2.10) that the initial cc pair (with transverse momentum l x in Fig. 1) is colorless. Indeed, for 
large values of k we can view the result of our calculation as a product of two factors: the probability to find a gluon 
(li) in the projectile hadron and its structure function in the target nucleus. The gluon structure function can be 
modeled by the interaction of a colorless probe such as dilaton or graviton ([33, 43]) with the nucleus through the 
splitting into the colorless cc pair. In Appendix we present a formal derivation of all the main results of this section 
by direct summation of the corresponding Feynman diagrams in the light-cone perturbation theory along the lines 
of the dipole model [33, 44]. 

Doing integrals over the longitudinal positions Zi of nucleons and summing over n in (2.10) we obtain 

r^R A 00 r 2R A r 2R A r 2R A 



/■2Ra 


/•2Ra r 


/ d Zi . . 


. / cfe n _i / 


'20 





= e Xp {-[a(x 2l r 2 ) + a(x 2 ,r' 2 )} P 2R A } ■ ^^,rS) P iR A _ ^ 
= exp {-a[x 2 ,{r- r') 2 ] P 2R A } - exp {-[a{x 2 , r 2 ) + a{x 2 , r' 2 )] P 2 R A ) . (2.11) 
Using (2.11) we can reduce (2.10) to the following expression 

da in (pA) _ XiG{x ^ m 2 c) [ d * rdze -ihkr f d 2 r > dz > yU-n' ' $ G (l u r,r' , z) 



dYd 2 kd 2 b 

x (exp{-a[x 2 ,(r-r / ) 2 }p2R A } - exp {-[a(x 2 , r 2 ) + a(x 2 , r 12 )] P 2 R A }) . (2.12) 



By forward rapidities we mean the direction of the projectile fragmentation. 
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This equation accounts only for the inelastic interaction and the physical meaning of (2.12) is the cross section of 
all possible inelastic interaction in which the cc pair is produced. We need to add the cross section for the elastic 
production of the quark-antiquark pair, which reads 

dvei(pA) =XiG{xum 2 c) f d 2 re ^k-r f d * r ' e^' $ G (h,r,r' , z = 1/2) 



dYd 2 kd 2 b 

x {l~cxp[-a(x 2 ,r 2 ) P 2R A }} • {l - exp[-a(a; 2 , r' 2 ) p2R A ]} . (2.13) 



The sum of (2.12) and (2.13) gives 



da tot (pA) _ XiG{xum 2 c) I d 2 re -ihh-t I d 2 r'e i ^'$ G (h,r,r , ,z = 1/2) 



dYd 2 kd 2 b 

x {1 - exp[-a(x 2 , r 2 ) P 2R A ] - exp[-<r(a;a, r l2 ) P 2R A ] + exp[-a(x 2 , {?- ?) 2 ) P 2R A }} . (2.14) 



Introducing the quark saturation scale Q 2 (see (A. 13) and (A. 14)) we can write 

a(x,r 2 ) P 2R A = 1 -r 2 Ql A (x). (2.15) 

The form of Q 2 is determined by the phenomenology of low x DIS [45, 46, 47, 48, 50] and forward hadron production 
in p(d)A collisions [49, 51, 52, 53]. Introducing a new dimcnsionlcss variable C = rn c r we can rewrite (2.14) as 

l^KiiOKtiO + KoiQKoiC) 
x {1 - exp (-C 2 Q'/4m 2 ) - cxp (~C' 2 Qll^l) + cxp [-(( - C? Q^m 2 ] } . (2.16) 

It was pointed out in Ref. [2] that the dominant contribution to the integrals on the r.h.s. of (2.16) is originating 
from the integration region r' <C r -C l/ni c) i.e. (' <C £ < 1 (or, equivalently, ( <C < 1). In this kinematic region 
(2.16) reduces to the following expression 

da tot (pA) 1 a s 2 
xiG(xi,m c ) 



dYd 2 kd 2 b An m 2 ir 



x dC 2 K (C)J (kC/2m c ) dC 2 K (C){l-exp[-C 2 Q 2 s (x 2 )/Am 2 c ]} (2.17) 
Jo Jo 

In the saturation region Q s 3> m c the dipole scattering amplitude reaches its unitarity limit 1 — e - ^ 2 s <( l2 '/ 4m c 1. 
Therefore, the rapidity distribution becomes 

da to t(pA) 2 . /OICA 

oc xiG(xi,m c ) ~ exp(— XY) , (2-18) 



dYd 2 kd 2 b 

while for the same process in hadron-hadron collisions we have (see (2.1)) 

datotipp) 



dYd 2 kd 2 b 



oc x\G(x\, m 2 ) X2G(x2, m 2 ) ~ constant(y) , (2-19) 



where we assumed that xG(x,m 2 ) oc l/x A at low x (which is true if Y is not too close to the proton fragmentation 
region). It is clear that there is a substantial difference between the rapidity distribution in these two cases. 
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Figure 3: The main diagrams for the process of inclusive cc production with fixed relative momentum in hadron-nucleus 
collisions (see Appendix). 



2.3 Nucleus- nucleus collisions in the KLN approach 



Nucleus-nucleus interaction can be characterized by the saturation scale which depends on the properties of both 
nuclei. In the KLN approach it is assumed that multi-particle production is entirely determined by the saturation 
scales of the colliding nuclei Q 2 S Ai (x\) and Q^a 2 ( x ^)- In the spirit of this approach we will generalize (2.14) to 
the case of nucleus-nucleus scattering using the Kovchegov's conjecture [54]. In [55] Kovchegov and Mueller noted 
that in order that their calculation of gluon production in pA collisions be self-consistent, an entire class of the final 
state interactions must cancel out in the light-cone gauge. Although they did not find a physical reason for that, 
Kovchegov suggested that the same conclusion may hold in AA collisions as well. Using this assumption he derived 
an expression for gluon production in heavy ion collisions in the light-cone gauge. In Appendix we derive (A. 33) 
along the same lines of reasoning. Here we would like to review the main steps. 

The main contribution to the inclusive cross section of cc production stems from the diagrams shown in Fig. 3. 
The sum of the diagrams of Fig. 3-a are proportional to 



(x-r') 2 



(2.20) 



while the diagrams of Fig. 3-b are proportional to 



xiG(xi,ml) 



l(*2) 



(2.21) 



The sum of (2.20) and (2.21) gives (2.14). 

For nucleus-nucleus collisions the main contribution stems from the set of the diagrams given in Fig. 4, which 




Figure 4: The main diagrams for the process of inclusive cc production with fixed relative momentum in nucleus-nucleus 
collisions (see Appendix). 



can be written as 

da tot (AA) _ 1 a s 4N C 



2 



dY d 2 k d?bd 2 V (2tt) 3 m 2 Tr a s ir 

d 2 (d 2 (' e <«-(f-?)/am. ^1 LL Kl{C)Kl{C) + KoiOKoiC)^ (2.22) 

i {l-exp[-C 2 QU/8m2]} {l-exp[-C 2 Q^ 2 / 8 m2]} 
+ ^ {l-exp[-C' 2 g 2 !yll /8m2]} {l-ex P K' 2 QU/ 8 ™c]} 

-JC^CT ^-^-^-O'QIaJ^I]} {l-exp[-(C-C') 2 Qs,A 2 / 8 ™c]}^ 

One can see that the first two terms in (2.22) are the same as (2.21) where factor xiG(xi,m 2 ) is replaced by (see 
(A.32)) 

^xM^mt) - (l < (2.23) 

or 

— xiG{x u m c ) ~^ ^ V 
while the last term in (2.22) is equal to (2.20) with the replacement 

a 7T 2 d 2 b ( (r-z') 2 Q 2 B 
^-x lG { Xl ,ml) — (l-e 
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R(Y) 
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Figure 5: The ratio R(y) defined in (2.27), for the process of inclusive cc production with fixed relative momentum k — 0) 
in nucleus-nucleus collisions. The calculation is performed for the Gold nuclei collision at RHIC energy ^/s = 200 GeV using 
the KLN value for the saturation scale [7]. 



To understand this replacement we notice that the last two lines in (2.10) can be written as 



2R A oo „2R A 

px 2 G{x 2 ,ml)dz e-^ x ^ + ^ r '^^ 2RA -^Yl / dz ^ 



^=0 J z a 



2R A *1R A 

dz n -i I dz n p n &i n (x2,r,r') 



2Ra 



px 2 G(x 2 ,m 2 ) dz exp < - <r(x 2 , r 2 ) + a(x 2 , r' 2 ) - ff m (i 2j r, r') p(2R A -z ) 



x 2 G(x 2 ,A/[r-r'] 2 ) 



a(x 2 ,r 2 ) + <j(x 2 ,r' 2 ) - a in (x 2 ,r, r') 



1 — cxp 



a(x 2l r 2 ) + a(x 2 y 2 ) - a in (x 2 ,r,r') 



p2R A 



a s TT 2 (r — r') 



A2 



1 — exp 



(r-r>) 2 Ql Ai (x 2 ) 



(2.25) 



(2.26) 



This corresponds to the diagram of Fig. 4-a. The sum in (2.25) reflects the fact that the dipole can scatter clastically 
only after the first inelastic interaction. Contribution of the diagram Fig. 4-b is treated in the same way. The low 
density limits, i.e. hadron-hadron or hadron-nucleus collisions are reproduced when Q 2 Ai <C m 2 and/or Q 2 A ^ <C m 2 . 
In Fig. 5 we plot the ratio 

da tot (AA) 



R(Y) 



dYd 2 kd 2 b 



k=0 



(2.27) 



da tot (AA) 
dYd 2 kd 2 b k=o,Y=o 

as function of rapidity Y. This ratio has a much sharper maximum at Y = than the corresponding ratio in pp 
collisions. 



3. J/tp production in hadron-nucleus collisions 
3.1 New production mechanism off nuclear targets 

The discussion and derivations of the previous section now allow us to turn to the main subject of our paper. Fig. 6 
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hadron — hadron collisions hadron — nucleus collisions 



Figure 6: The process of inclusive J/ij) production in hadron-hadron (Fig. 6-A) and in hadron-nucleus collisions (Fig. 6-B). 

displays the J/if> meson production in pp and pA collisions at the leading order in a 2 A 1 / 3 . The cross section for the 
latter is a direct generalization of (2.1) and reads 3 

2R A z 

dcrjpA) „ f f d% 

dYd 2 b 



<x / pdz Q / pdz x / _l0 G ( a . lj | 1 )2 d 2 rdz^ G {l 1 ,r,z)®^ v {r,z) (l - e^) (l - e^) 



o o 



x2 f d 2 r'dz'$* G (h,r',z)®iZ>* v (r',z') (l - e'^^ (l - e'^ 3 ^ 



^MX2,12) J ^<t>G(X2,k) (3.1) 

where ® VEV is projection of the J/ip light-cone "wave-function" onto the virtual gluon one. Trace over all 
relevant quantum numbers is implied. Assuming that ^V(r, x) oc S(z — 1/2) S(r) and I 2 z(l — z) <C m 2 , i = 1,2, 3, 
this projection takes the following form [2, 32] 



'3rj/,,, . B+e - M.r • 



vMm c ,r,z) ®*v(r,*) = l f~ J/ tr +< ^ "f^Kr) • (3.2) 

y 48 7r a em 4 

At short distances r, r' < l/m c the cross section in (3.1) is proportional to r 2 r' 2 . This fact reflects the higher 
twist nature of the suggested mechanism. Since the position of the pair of nucleons is not fixed the full contribution 
should be proportional to A 2 ! 3 while the mechanism of Fig. 6-A leads to an enhancement by A 1 / 3 . The enhancement 
factor stems from integrations over z and Z\ in (3.1) 

[ RA p^[° P ^ = \ {p 2R A f. (3.3) 

Parametrically, the mechanism in Fig. 6-B is different from that in Fig. 6-A by the factor agA 1 ^ 3 . Therefore, in the 
spirit of the quasi-classical approximation in which we assume that a 2 A 1 / 3 ~ 1 we conclude that the mechanism in 
Fig. 6-B is enhanced by a big factor l/a s . 

3.2 Propagation of the colourless cc pair through a nuclear target 

In Ref. [2] detailed arguments were given which justify the application of the dipolc model for calculation of J/ip 
production at forward rapidities at RHIC. It has been argued that the coherence length for the cc pair is sufficiently 

3 All cross sections in this and the following sections are for the J/ip production. 
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larger than the longitudinal extent of the interaction region. This means that the development of the light-cone 
"wave function" happens a long time before the collision. We relied on this physical picture in our calculations in the 
previous section. Concerning J/ip, its production is characterized by an additional formation time scale proportional 
to the inverse binding energy ~ a 2 m c ; for a quantitative estimate, see [56]. This time is certainly larger than the 
charm quark pair production time (by a factor of about 1/a 2 ) implying that the formation process takes place far 
away from the nucleus. Therefore, in the following, we concentrate on the dynamics of cc pair interactions with the 
nucleus. 




Figure 7: The process of inclusive J/ip production in hadron-nucleus collisions due to the interaction with an even number 
of nucleons. 

Since the soft gluon emission processes are suppressed in the quasi-classical approximation, the J/ip meson is 

predominantly produced through the hadronization of the color singlet cc pair. This rules out diagrams of the type 

Fig. 3-b corresponding to the clastic interaction. The diagrams of the type Fig. 3-a are shown in Fig. 7. Note that 

since the J/ip quantum numbers are 1 while those of gluons are 1 _ an odd number of gluons must connect to the 

charm quark line. Consequently, each inelastic interaction of the cc pair must involve two nucleons. To take this into 

account we write an analogue of (2.10) in which the sum over all inelastic processes (i.e. sum over n) involves only 

even number of interactions. We have 

da in (pA) 2 
dYd 2 b = C F xiG{x u m c ) (3.4) 

x / pa in (x 2 ,r,r')dz / d 2 r * G (h, r, z = 1/2) * v (r) <g> / dV^ft,r', Z = l/2)^(r') 



o 



w r 2R A f 2R A r 2R A 

<^) + a( x ^)) P 2R A J- dZl dz 2 ... dz 2n+lP 2n+1 a^( X2 ,r,r') 

,v,— n J Zi) J Z\ J Zn„ 



x I e 

The color factor in front comes from the calculation of Fig. 6-B, namely, it is equal to 

Tv{t a t b t c ) Tv{t a t b t c ') 5 CC ' - 4 S ° C ' f-"c + dabc <W ) (3.5) 

16 

(N 2 - 1) ( N N 2 ~A \ = N 2 -l N 2 -2 = / iV r 2 -l \ 2 N C (N? - 2) ^ ^ 



16iV c y c ' N c J 2N C 4 \ 2N C J 2{N 2 - 1) ' 

Since <7j„ is proportional to Cf = (N 2 — l)/2N c , we extract this factor from the color coefficient of (3.5). The last 
of equations in (3.5) is written in the large N c approximation. 
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We argued in Sec. 2.3 (see Fig. 3) and in the Appendix that cc pair in the color octet state passes through the 
target with the same clastic (sec Fig. 3-b and (2.21)) and inelastic (see Fig. 3-a and (2.20)) cross sections as the cc 
pair in the color singlet state. This is the reason we do not need to change (3.4) to include the color octet state 
interaction with the target. 

After integration over Zj's and summation over n using the identity Yl^Lo a2 ™ +2 /(2n + 2)! = cosh a — 1 we obtain 
the following formula 

da in (pA) 
dYcPb 



= C F x l G{x u m 2 c ) / d 2 r* G (h,r,z = 1/2) tfy(r) / d 2 r'* G (h, r',z = 1/2) ^* v {r') 

exp - (a(x 2 ,r) + a(x 2 ,r') + a in {x 2 , r, r')) p2R A 
2 exp - (cr(x 2 , r) + cr(x 2 , r')) p2R A 



x 2^ exp 



v(x2,(r-r') 2 ) P 2R A 



(3.6) 



The color factor in (3.4) as well as in (3.6) corresponds to the diagram of Fig. 6-B. 

3.3 J/\& production in hadron- nucleus collisions in the saturation regime 

In the quasi-classical approximation the gluon saturation scale is given by [33, 43], see (A. 13) 

Ql A {x) = A7T 2 a 2 sP T(b), 



(3.7) 



where p is the nucleon density in a nucleus, N c is the number of colours, b is the impact parameter and Tib) is the 
optical width of the nucleus. Eq. (3.7) determines the scale of the typical transverse momenta for the inclusive gluon 
production [44]. Its value was extracted from the fit to the hadron multiplicities in nuclear collisions at RHIC [10, 7]. 
However, for the penetration of the quark-antiquark pair the typical saturation scale is about twice as small and 
we refer to it as the quark saturation scale Q 2 A , see (A. 14). This scale was extracted from fits of the F 2 structure 
function in DIS [43, 45, 46, 47, 48] as we have already mentioned. Both phcnomcnological approaches agree with 
each other, so the use of either quark or a properly rescalcd gluon saturation scale is merely a matter of convention. 
In this paper we will use the gluon saturation scale (3.7). Using this definition for the saturation momentum we have 
o(x 2 ,r) = r Q 2 (x 2 ) /8 (cf. (2.15)). Substituting this expression into (3.6) we can re-write it in a more convenient 
form 



da in (pA) 
dYd?b 



Cp xiG{xi 1 m 2 c ) 

x { d 2 r^ G {h,r,z= \/2)®^ v {r) [ d 2 r' V G (h, r', z = 1/2) ® $* v (r') 



(3.8) 



x - exp 



exp 



2 exp 



{r 2 + r l2 )Q 2 s 



Integrating over the angle between r and r' we derive 
da in {pA) _ N c (N 2 -2) 



dYd 2 b 



2(N 2 -1) 



x\G(xi, m c ) 

2 r^ G (h,r,z = 1/2) ^ v (r) ® / dV * G (h, r',z = 1/2) ** v (r') 

h I 



(3.9) 



x exp 
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Deeply in the saturation region where Q s .a 3> m c the typical dipole sizes are much smaller than l/m c . Thus, 
we can expand the wave function (3.2) ® ~ const. The main contribution comes from (r — r*) 2 < l/Q 2 ^ 
while r \jm c . It gives 

do-<„(pA) oj^cr^^a^Q^^j a cxp (_2AF) . (3.10) 



dFd 2 6 



In deriving (3.10) we used the same assumptions as in the case of cc-pair production with fixed relative momentum 
(see (2.18)). One can see that (3.10) leads to a rapidity distribution that is more narrow than the distribution in 
hadron-hadron collisions given by (2.19). 



4. Inclusive J/ip production in nucleus— nucleus collisions 



Using the same arguments as in Sec. 2.3 which led us to (2.22) we can generalize (3.9) to obtain our main result 
the formula for J ftp production in nucleus-nucleus collisions. It reads 

1 dtr(AA) C 2 F 



S A dYd 2 b 4tt 2 c 



x 
2r 



d 2 r % G {h,r, z = 1/2) ® # v (r) / d 2 r' % G {h,r', z = 1/2) ® %* v {r') (4.1) 
l — jcxp (-gfc- r') 2 (Ql Al + Q^ 2 )) - CX P (- \( r -+ r -') 2 (OU + QIaS) 

- exp (-j(^r') 2 QU ~ \{r 2 + r' 2 )Ql A )j + exp (^r + /) 2 QU - \{r 2 + r' 2 ) 

- exp (4(r-r') 2 QU - I(, 2 + / 2 )Q 2 ^ + exp (-Ife + r') 2 Q 2 A2 - I(r 2 + r' 2 ) Q 2 ,^) } , 

where is the transverse overlap area. One can check that this formula describes the hadron-nuclcus J / ip assuming 
that Q 2 Ax ^ s small- 




Figure 8: The process of inclusive J/ip production in nucleus-nucleus collisions due to inelastic interaction with both nuclei. 

In deriving (4.1) we summed up the inelastic cross sections for both nuclei. Let us denote the nucleon coordinates 
in the nucleus A\ by z\ and in the nucleus A-^ by z\. The number of possible inelastic interactions (see Fig. 8) is odd 
for both nuclei. For a fixed number of total inelastic interactions In — 1 the number of interactions in each nucleus 
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can be 1 < k < 2n — 2. Therefore, for the term inside the curly brackets in (4.1) we have 



l ^ So M {l Q2s ' A2 ) (iQlA.) (^■r') 2 dz a dz' exp^- 1 -(r 2 + r' 2 )(Ql Ai +Ql 

2n-2 ,2Ra 2 p2Ra 2 ^Ra 2 /j \ k-1 

x dzi dz 2 ... dz k _ lP k [^QIa 2 2 L- 

k=l z ° Zl J z k-i V 

°° t 2R Al r2R Al r2R Al /-, \ 2n-fc-2 

E /, dz l , dZ 'f - , ^n-k-lP 2 "-* -1 U '^ 2 -'- ^ 



X 

n=2 " z o 

Using the following mathematical identity 



with j = 27i—l wc obtain 



l /i \ 2n— 1 1/1 \ 2n— 1 



n=2 



-, v 2?i— 1 ' 



(2n-l)! V8 

= sinh Q (Q» Al + Q2 A J 2r • r'J - sinh Q Q* Al 2r ■ r') - sinh Q Q 2 .^ 2r • r') . (4.4) 

yielding (4.1). If Q 2 S A<2 <C Q 2 A we can expand (4.4) using sinh(a + b) — sinh a — sinh b ss b (cosh a — 1) + 0(& 2 ); then 
(4.1) reduces to (3.8). 

For a qualitative discussion it is instructive to rewrite (4.1) in the region r' <C r ~ l/Q s ,A "C l/m c . Expanding 
expression in the curly brackets we derive 

{•••} = ~ QU QIa 2 (Qi Al + Ql,A 2 ) (r ■ r_r (1 + 0(Ql A r' 2 , Q% A r' 2 r 2 )) cxp (-~(Q 2 , Al + Q% M ) r 2 ) . (4.5) 

In this approximation Eq. (4.1) becomes (after integration over the angle between r and r') 
1 da(AA) 



; cx / d i r* G {h,r,z= 1/2) ®* v {r) J d 2 r' ** G (h, r' , z' = 1/2) ® ** v (r') 

*Ql Al Ql,A 2 («U + Qs,A 2 ) ^ r a exp {-r 2 (Q 2 Ai + Q^J/S} (4.6) 



(4.7) 



In (4.7) we replaced the wave functions of (3.2) by a constant as already discussed in Sec. 3.3 and took the integral 
over the angle between r and f . 

From (4.7) one can see that the spectrum of J/^'s in ion-ion collisions is more narrow than the one in hadron- 
hadron. Explicitly 

da(AA) da( PP ) 1 dajpp) 3MY \ , Af * 

dY * d Y (Qi Al + Qi A2 r dY ' (4 - 8j 
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Figure 9: Rapidity dependence of the ratio R(Y) 
the KLN expression was used. 



#(y=0) 



for the gold-gold collision at RHIC. For the saturation momenta 



where we use that Q 2 (x) oc and Y is the rapidity of J/ip in the center-of-mass frame. Therefore, our prediction 

is that the rapidity distribution of J/ip is much more narrow in nucleus-nucleus collisions than in the proton-proton 
ones. 

To evaluate how close we are to the saturation region at RHIC energies in this process we first rewrite the general 
formula for the kinematic region r 3> r' . It takes the following form: 



1 da (A A) 
~Sa dYd 2 b 



^QiAM)Qi,A^){Qs,AM) + QUA^)] / <KCK 2 (c)e 



[QI Ai ^i) + QIa 2 ^2)} 



(4.9) 



In Fig. 9 we plot the result of our calculation for the ratio R(Y) = jy(Y)/jy(Y = 0) using (4.9). For Au - Au 
collision at RHIC with y/s = 200 GeV taking the KLN value for the saturation momentum Q 2 (y = 0) = 2.2 GeV 2 
for central collisions we find that the rapidity distribution turns out to be very narrow although not quite to an 
extent suggested by the approximate expression (4.8) (see Fig. 9). The rapidity distribution in Fig. 9 is driven by 
the ratio Q 2 Ai Q 2 A2 /{Q 2 a 1 + Q 2 a-,) 3 - The cross section decreases with the increase of the value of the saturation 
momentum Q 2 Ai + Q 2 A . However, this decrease is much milder than in (4.8). 



5. Numerical calculations 



In this section we perform numerical calculations of inclusive J/ip production using (4.1). First of all. we reinstall the 
impact parameter dependence of the saturation scales and consider a realistic distribution density for nuclei. Recall 
that Q 2 oc pT(b). Denote the impact parameter between centers of two nuclei as b. The position of a nuclcon inside 
nucleus A\ with respect to its center denote by s. Then the position of a nuclcon in the nucleus A 2 is given by b — s. 
We have 



ls,A 



Qs,A 2 



Qs,As 



(k-s) 



(5.1) 



In our Glauber-type approximation (see e.g. [57]) we neglect the impact parameter dependence in nucleon-nucleon 
interactions considering their range much smaller than the size of nuclei. The observable that we are going to 
calculate is the number of J/V' s inclusively produced in nucleus nucleus collisions at a given rapidity Y and a 



- 15 - 



350 



300 - 




-6 -4 -2 2 4 6 

y 



Figure 10: J ftp rapidity distribution in Au-Au collisions for different centrality cuts. Experimental data from [59]. 



ccntrality characterized by the impact parameter b. The corresponding expression reads 

dN AA (Y,b) 
dY 



oc / d 2 s Q 2 sAl (an , s) Q 2 sM (x 2 ,b- s) [Q 2 StAl (an , s) + Q 2 sM (x 2 ,b- s)] 

J d(C 9 K 2 (() ex P (-^[Ql Al (x u s) + Ql A2 (x 2 ,k- s)}} . (5.2) 



i ai\ (X-0 aa (1,0) n 

where = dy '' and 



xi = — , x 2 = — (5.3) 
V s V s 



We can also write down (5.2) in the following way 
dN AA (Y,b) dNPP(Y) 



°~ dY 



[ d 2 s T Al (s)T A2 (b- s) [Ql a {xi,s) + Q 2 A {x 2 ,b- s)] -X, 



J d( C 9 K 2 (C) exp | - ^ [Ql Al (an , s) + Q% M (x 2 , b - s)} | 



(5.4) 



The overall normalization constant C includes the color and the geometric factors Cp/(4Tr 2 a s S p ) where S p is inter- 
action area in proton-proton collisions; it also includes a rather poorly known amplitude of charm quark-antiquark 
transition into J/ip and a gluon in the case of pp collisions (see Fig. 6-A). 

To calculate the multiplicity of J/ip's in AA collisions using (5.4) we need to know (i) rapidity distribution of J/ip 
multiplicity dN pp /dY in pp collisions and (ii) the overall normalization constant C. We fitted the rapidity distribution 
of J ftp's in pp collisions to the experimental data of Ref. [58] with a single gaussian. The global normalization factor 
C is found from the overall fit. 

Now we can compute the nuclear modification of the J fip rapidity distribution in ^4^4 collisions at all centralitics 
using (5.4). In figure 10 we compare our results with the experimental data of PHENIX Collaboration [59] for Au-Au 
collisions at ^/s = 200 GeV. The agreement between our calculation and the experimental data is reasonable. 

In figure 11 we show the same results in the "measured/expected" form, namely the experimental data divided by 
our calculations. It is tempting to use the difference between this ratio and the unity as a measure of the magnitude 
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results. 



Ratios of J /if) rapidity distribution in Au — Au collisions: Experimental data [59] divided by theoretical our 



of the final state effect for J/tp production. However, at the moment we prefer to refrain from a premature conclusion 
about the size of the final state effect as such a conclusion would crucially depend on the value of the C factor. In 
addition we are aware of the limitations due to the accuracy of the experimental data. To give a firm conclusion we 
need to perform a comparison with a high precision dAu data, where the C factor can eventually be fixed with a 
higher accuracy. 

To emphasize the nuclear dependence of the inclusive cross sections it is convenient to introduce the nuclear 
modification factor 

dN AA 

Raa(v, N part ) = dv dNPP . (5.5) 

^VCOII -ay- 
it is normalized in such a way that no nuclear effect would correspond to Raa = 1- In Fig. 12 we plot the result of 
our calculation. The nuclear modification factor exhibits the following two important features: (i) unlike the open 
charm production, J/ip is suppressed at y = 0. This is not very surprising since the probability of the J/ip formation 
is reduced due to multiple interactions with the gluons; (ii) cold nuclear effects account for a significant part of the 
J/tp suppression observed in heavy ion collisions. On the other hand, the cold nuclear effects discussed in this paper 
may not be sufficient to account for a very low Raa m the most central events. Higher precision dA and AA data 
will allow to tell whether there is a suppression of J /if) in quark-gluon plasma, or the directly produced J/ip's survive 
[60]. 

6. Conclusions 

In this paper we have developed a model for J/-0 production in heavy ion collisions. Assuming that contributions 
of strong color fields of the two colliding nuclei do not interfere we summed all the diagrams proportional to the 
positive powers of the large parameter a^A 1 ^ 3 ~ 1 in both nuclei. Our main result is given by Eq. (4.1). In the 
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N Da rt 

Figure 12: Nuclear modification factor for J/ip production in heavy ion collisions for different rapidities. 

RHIC kinematic region this formula can be simplified and is given by (4.9). We used this equation in our numerical 
calculation with the realistic nuclear profiles described in Sec 5 in details. The results are presented in Figs. 10-12. 

We can see that the rapidity and centrality dependence are reproduced quite well. This observation implies 
that an appreciable amount of the J/ip suppression in high energy heavy ion collisions comes from the cold nuclear 
effects. Fig. 12 demonstrates that the nuclear modification factor for J/ip production is strongly suppressed even at 
zero temperature. While J/ip suppression in the forward direction is not a surprise (the nuclear modification factor 
for light and, perhaps, heavy hadrons is known to be suppressed), similar behavior in the central rapidity region is a 
peculiar feature of the J/ip production. The reason is that multiple scattering of cc pair in the cold nuclear medium 
increases the relative momentum between the quark and antiquark, which makes the bound state formation less 
probable. Formally, the sum rule proven in [11], which guarantees emergence of the Cronin enhancement for single 
partons, is broken for the bound states. 

Our result strongly suggests that the cold nuclear matter effects play a very important role in J/ip production in 
heavy ion collisions. The final nuclear modification factor, which is measured in experiment, is undoubtedly a result 
of a delicate interplay between the cold and hot nuclear matter effects. 

We consider the present work as the first step towards understanding the role of the cold nuclear effects in 
J ftp production in high energy heavy ion collisions. We calculated the parametrically enhanced contribution coming 
from even number of scatterings of cc pair in the nuclei. However, we neglected other contributions that may be 
phcnomcnologically important though parametrically small. These include soft gluon radiation in the final state and 
color octet mechanism of J/ip production. Moreover, for peripheral collisions these contributions become of the same 
order as the one discussed in this paper. Therefore, we plan to perform a detailed investigation of these contributions 
in the future. 
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A. Heavy quark production 




Figure 13: Diagrams contributing to the qq production in pA collisions at high energies. Transverse momenta of the produced 
quark and antiquark k lt k 2 , their transverse coordinates in the amplitude x 1 , x 2 and in the complex conjugated one y 2 
and the gluon transverse coordinate in the amplitude u and in the complex conjugated one v are shown in boldface. Vertical 
dashed lines indicate all possible interaction times of incoming proton with the nucleus (shown by crosses). More detailed 
discussion can be found in [36, 37]. Notations follow [37]. 



In this appendix we derive the cross section for production of a cc pair with fixed relative momentum in pA 
collisions. A general problem of qq pair production in pA collisions including all possible nonlinear evolution effects 
was solved in Refs. [36, 37]. Let us introduce the following notations: k t and k 2 are the produced quark and 
anti-quark transverse momenta, q is the gluon transverse momentum, z = k + /q+ is a fraction of the light-cone 
momentum of gluon carried by the produced quark; £j and y are the transverse coordinates of the produced quark 
in the amplitude and in the complex conjugated amplitude respectively; x 2 and y 2 are the corresponding coordinates 
of the antiquark. Transverse coordinates of gluon in the amplitude u and in the complex conjugated amplitude v are 
given by u = zij + (1 — z)x 2 with u = \u\ = r and x 12 = X\ — x 2 = Tj (#12 = |x 12 |) and analogously for v. With 
these notations the double inclusive quark-anti-quark production cross section is given by [36, 37] 



da 



d 2 k 1 dk 2 dyd 2 b 4(2tt) £ 



d Xi d x 2 d yi d y 2 j dz e 




-i k 1 -(x 1 -y i )-ik 2 -(x 2 -y 2 ) 



(A.l) 



x y ] (—1 ' —2 5 Vjy ' —2 5 ^0 '—'ij (—1 ) —2 ' Vjy ! —2 ! ^) 

The products of the light-cone "wave functions" are detailed as follows [37] 

2 

^nfe, x 2 ;y v y 2 ,z) = AC F 



F 2 (x l7 x 2 ; z) F 2 {y y ; z) [(1 - 2 z)< 

- 1 - 2 xi 2 yi2uv 



(A.2) 



x (x 12 ■ u) (y i2 ■ v) + (dj UiX 12 j) (e k i v k y 12 i)} + F 1 i(x 1 ,x 2 ; z) F 1 {y l ,y^, z) m 2 



+4z 2 (1 - zf Fofe,^; z) F (y v y 2 , z) - 2z (1 - z) (1 - 2 z) 



X\2 U 



F 2 (x 1 ,x 2 ;z) 



x Fo(y v y 2 , z ) 



y ■ v 

2-12 - 

yi2v 



F2(y v y 2 ,z)F (x 1 ,x 2 ;z) 



(A.3) 
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$22 fe,, x 2 ;y y ;z) = AC F ( — J to 2 (^(minjiCifmyn) , 

- 1 - J \ 7T / ^ ^12 2/12 m v 

x (x 12 • u) (y 12 • + (eij UiXuj) (efci Wfc yi2;)] + ^o(m^i2) K (my 12 ) 



[(l-2z) : 



*i2fe,^2;y i; y 2 ; 2 ) = -4CV I — I F 2 {x 11 x 2 ;z)Ki(myi 2 ) 



1 



•^12 2/12 



[(l-2z) ; 



(A.4) 



x (^12 -m) (y 12 ■«) + (eijUiX 12] ) (e k i v k 2/i 2 ;)] + m F^^, x 2 ; z) K (my 12 ) 

-2 z (1 - *) (1 - 2 z) Fofe,^; z) #i (m Wa ) 1 , 

yi2 ^ J 

*33fe,^2;y i; y 2 ; z ) = ^life)^;^)^;*) + $22(^1, ^y^; 2 ) + ^^fe,^;^)^'^ 

+$2i {x^x^^yj z) 
^nix^x^y^y^z) = -^iiix^x^y^y^, z) - $i 2 fe,x 2 ; y v y 2 ; z) 
^2z{x l ,x 2 ;y v y 2 ,z) = -^2i(x 1 ,x 2 ;y_ 1 ,y_ 2 ; z) - ^ 22 (x 1 ,x 2 ;y_ 1 ,y 2 ; z) 
®i J \x 1 ,x 2 ,y v y 2 \z) = ^{VvVyX^x^z). 
The auxiliary functions F\ , F 2 and Fq are defined as 



F 2 (x 1 ,x 2 ;z 
Fx^x^z) 
Fq{x 1} x 2 ;z) 



) = J dqj 1 (qu) Ki ^x 12 \J m 2 + q 2 z(l — z)j \J m 2 + q 2 z(l — z) , 



dq Ji (q u) K \X\2 \Jm 2 + q 2 z(l — zjj , 
dq q J (q u) K ^x 12 ^ m 2 + q 2 z(\ - zj^j 



where u — x 12 = x x — x 2 , x\ 2 = \x 12 \, and q = k 1 + k 2 . 
By definition, the gluon saturation scale Q s 

Ql = 47ra 2 s pT(b) 



(A.5) 

(A.6) 
(A.7) 
(A.8) 
(A.9) 

(A.10) 
(A.ll) 
(A.12) 

(A.13) 



with p the nucleon number density in the nucleus and T(b) the nuclear profile function. We also use the quark 
saturation scale Q 2 give by 

Ql= ^-Ql^'-Ql (A.14) 

Throughout the theoretical discussion we assumed for simplicity that the nuclear profile is cylindrical with T(b) 2Ra 
for b 2 < R 2 A . However, numerical calculations in Sec. 5 are performed with an accurate parameterization as discussed 
there in detail. In the large iV c approximation we write 

(A.15) 



^ii(x_i, x^y^y^z 
^22{x_i, x 2 ,y v y 2 ,z 
^■i{xi, x 2 ;y v y 2 ,z 
^■i2{x 1 ,x 2 ;y v y 2 ,z 
^Ax\iX 2 \y_ v y_ 2 \z 
^is{xi, x 2 ;y v y 2 ,z 



e -| (xi-V^ 2 Ql InCVI^-yJ (x 2 -y 2 ) 2 Q 2 t \n(l/\x 2 -y 2 \ 

e -\ (m-v) 2 Q 2 B \n(l/\u-v\ m) 



e -| (Hi-V.) 2 Ql ln(l/| £l -3i| {x 2 -h) 2 Ql ln(l/|x 2 - £ | n) 

e -\u 2 Ql ln(l/tt M ) : 

„-5i?Q' lnCl/xt^-i^ ln(l/x 2 /i) 



(A.16) 
(A. 17) 
(A.18) 
(A.19) 
(A.20) 
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All other Cj/s can be found from the components listed in (A.15)-(A.20) using 

^jfe.^;^,^^) = ^ 3l (y v y 2 ,x 1 ,x 2 ;z) (A.21) 

similar to (A. 9). 

If the typical gluon momentum q is much smaller than the produced quark mass, the above expressions can be 
significantly simplified. Indeed, since z(l — z) < 1/4 we get 

q 2 z(l - z) < m 2 , (A.22) 

and the auxiliary functions read 

F2(x 1 ,x 2 ;z) — Ki(x\2m) m u _1 , (A. 23) 

Fi(x.i,x 2 ; z) = K Q (xi 2 m) u^ 1 , (A.24) 

F (^,x 2 ;z) =0. (A.25) 

In this approximation the only non-vanishing products of "wave functions" are given by 



2 



4C F — m 2 \ K^m) K^m) 2 -^[(l-2zY(x 12 -u)(y.-v) 

7T / L a;i2yi2M 2 « 2 Ji 



li ■ V 

+(eijUiX 12 j) (e k i v k yi2 1)} + Xo(a;i2m) K (y 12 m) =^ 



Averaging over all directions of gluon emission from the valence quark using (ey u,; X12 j tki Vk y\2i) = uwi^ • y 12 we 
arrive at the well-known result [36, 40] 

$ 11^, x 2 ,y v y 2 ,z) = 



4C F — — <^ ^^[{l-zY + z 2 ]K l {x l2 m)K 1 {y l2 m)+K {x l2 m)K {y l2 m) \ (A.26) 
Let's now introduce the following notations, see Fig. 1: 

p=k 1 + k 2 , k = k x - k 2 + (1 - 2z) (k x + k 2 ) , (A. 27) 

or, cquivalcntly, 

h=zp+h, k 2 = (l-z)p-~k. (A.28) 

At z = 1/2, momentum k becomes the relative momentum of the cc pair. In this case integration over the total 
momentum p results in the delta function 8^ 2 \u — v). This delta-function simplifies expressions in the exponents of 
(A.15)-(A.20). For the sum over all rescattering factors including the signs of Qij's we get (omitting the logarithms 
ln(l/x/z) for brevity) 

2fe 2 ,y 12 ; z=l/2) = e -5(iEi2-M 12 ) 2 Q" + 1 _ e -h£ a Ql _ e -\Q 2 ,vl 2 (A.29) 

Using the delta function in (A. 2) to integrate over u and integrating over v in the leading logarithmic approxi- 
mation we derive 



da CFCtlm 2 



d 2 kdyd 2 b 4tt 5 



dz / d 2 x 12 / dW-'H-feu-Hia) ln(l/ M i £l2 - yi J) 
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~ 12 [(1 - zf + z 2 ] Ki(x 12 m) Kx(y l2 m) + K Q (x 12 m)K (y 12 m) 1 3( £l2 ,y ; z) (A.30) 
xi 2 y\ 2 J - 12 

Finally, using the x 12 = r, y± 2 = r' notation obtain the final expression for the cc pair production with fixed relative 
momentum (at z — 1/2): 

d2 d k a d P y A d2b = ^ J d 2 r J d 2 r' e-^U-fe-r') x G( Xl ,mt){~^ K^rm) K^r'm) + K (rm)K (r'm) 

N < e-lfe-r') 2 Q 2 S + ! _ e -fr 2 Q 2 _ e -|Q>' 2 | (A.31) 

To generalize this expression in the case of heavy ion collisions we replace the proton's gluon distribution function 
by the correlation function of the color field potentials of the second nucleus as follows [54] 



d 2 b ~K 2 a s r 2 

2 ~ ir>2 



l-e-i rQ A. (A.32) 



Therefore, taking the large N c approximation Q s f=a and GV « N c /2, we derive 
diV AA « s m 2 f 2 f 2 -4fe.(r-r') f 1 E ■ r' 



(Tr / d 2 r' e- l ^ (r -~ L H-^^ K 1 {rm)K 1 (r'm) + K (rm)K (r'm) 



d 2 kdyd 2 b Sir 4 J J [2 rr 

8C F f 1 



(l _ e-^ 2 *) (l - e-^Q^ + ^ (l - < -' ^' ) (l - e-fe' 2 «-) 



(r — r') 2 



1 — e s 



) (i- e -ife-r') 2 Q? 2 ) j. (A.33) 
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